*****************************************************************************
* Code for generating weather data at the county-year level
*
* - Based on file "degreeDays.do" from Schlenker's website
* - Modified to generate counterfactual temperatures
* 
*****************************************************************************

use "Data_new/gridInfo.dta", replace
levelsof(fips), local(fipsList)

******************************************************************************
* loop over fips
local appendID = 0
foreach s of local fipsList {
    * loop over years
	forvalues t = $yearMin/$yearMax {
        ******************************************************************************
		display("Working county `s' in year `t'")
        * load data
    	use if (dateNum >= mdy($monthBegin,$dayMonthBegin,year(dateNum))) ///
    			& (dateNum <= mdy($monthEnd,$dayMonthEnd,year(dateNum))) ///
    			using "$path_rawData/year`t'/fips`s'.dta", clear
        
		* Replace min and max temperature adding value under the counterfactual (0 if baseline)
		replace tMin = tMin + $counterC
		replace tMax = tMax + $counterC
		
		gen tAvg = (tMin+tMax)/2
        label var tAvg "average temperature - average of minimum and maximum temperature"

        * merge in grid weights
        * make sure there is no (_merge == 1) as gridInfo should have data for each gridNumber
        qui merge m:1 gridNumber using "Data_new/gridInfo.dta", assert(2 3)

        * keep only were weights are nonzero
        qui keep if (_merge == 3) & (cropArea > 0)
        drop _merge

        if (_N > 0) {
            ******************************************************************************
            * generate degree days
            foreach b of global boundList {
                * create label for negative bounds by adding Minus
                if (`b' < 0) { 
                	local b2 = abs(`b') 
                	local bLabel Minus`b2' 
                } 
				else { 
                	local bLabel `b' 
                }

                * default case 1: tMax <= bound
                qui gen dday`bLabel'C = 0

                * case 2: bound <= tMin
                qui replace dday`bLabel'C = tAvg - `b' if (`b' <= tMin)

                * case 3: tMin < bound < tMax
                qui gen tempSave = acos( (2*`b'-tMax-tMin)/(tMax-tMin) )
                qui replace dday`bLabel'C = ( (tAvg-`b')*tempSave + (tMax-tMin)*sin(tempSave)/2 )/_pi if ( (tMin < `b') & (`b' < tMax) )
                drop tempSave
            }

            ******************************************************************************
			* save variable labels
			foreach v in tMin tMax tAvg prec {
				local l_`v': variable label `v'
			}
             
            * collapse by aggregation level
            collapse tMin tMax tAvg prec dday* [aweight=cropArea], by(fips dateNum)

            * collapse by year
            gen year = year(dateNum)
            collapse (sum) prec dday* (mean) tMin tMax tAvg, by(fips year)

			* label variables
			foreach v in tMin tMax tAvg prec {
				label var `v' "`l_`v''"
			}
			foreach b of global boundList {
				* create label for negative bounds by adding Minus
				if (`b' < 0) { 
					local b2 = abs(`b') 
					local bLabel Minus`b2' 
				} 
				else { 
					local bLabel `b' 
				}
				label var dday`bLabel'C "total degree days above `b'C (Celcius and days)"
			}
			label var year "year"

            ******************************************************************************
            * save data
            if (`appendID' > 0) {
	            append using $save_name
	        }
	        local appendID = 1
            sort fips year
            qui save $save_name, replace
        }
    }
}

* clear up files
! rm gridInfo.dta
